setwd("/Users/chuangchen/Library/CloudStorage/OneDrive-UniversityofPittsburgh/Pela project")
library(haven)
library(dplyr)
library(ggplot2)
library(ggpubr)
options(scipen=999)


overall <- data.frame(left=c(0.2626, 0.2924), right=c(0.1170, 0.1109))


for(i in 1:2){
  data_name <- paste0("model", i, "_pred_prob.dta")
  
  data1 <- read_dta(data_name) %>% filter(`_predict`==1) %>%
    select('_at1', '_margin', '_ci_lb', '_ci_ub')
  E <- c("Ec","Em","Ep")
  data1$eco <- rep(E, 3)
  M <- c("Mc","Mm","Mp")
  data1$moral <- rep(M, each=3)
  colnames(data1) <- c("eco_moral", "pred", "ll", "ul", "eco", "moral")
  
  left <- ggplot(data1, aes(x=eco, color = moral)) + 
    geom_vline(aes(xintercept="Ec"), linewidth=30, color = alpha("gray", 0.2)) +
    geom_vline(aes(xintercept="Em"), linewidth=30, color = alpha("gray", 0.2)) +
    geom_vline(aes(xintercept="Ep"), linewidth=30, color = alpha("gray", 0.2)) +
    geom_pointrange(aes(y=pred, ymin=ll, ymax=ul), position=position_dodge(0.5)) +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"),
          axis.text = element_text(size = 12),
          axis.title.y = element_text(size = 14),
          axis.title.x = element_text(size = 14),
          axis.text.x = element_text(color = "black"),
          axis.text.y = element_text(color = "black")) +
    geom_hline(yintercept=overall[i,1], linetype="dashed", color = "red") +
    annotate("text", x=1, y=overall[i,1]+0.02, label=paste0("overall probability=", round(overall[i,1],2)), size=3, color="red") +
    ylab("Predicted Probability of Being Leftist") +
    xlab("Economic Position") +
    ggtitle("Left") +
    scale_x_discrete(breaks=c("Ec", "Em", "Ep"),
                     labels=c("conservative (Ec)", "middle (Em)", "progressive (Ep)")) +
    scale_y_continuous(limits = c(-0.1,0.9), expand = c(0, 0)) +
    scale_color_manual(values = c("Mc"="#FF99FF", "Mm"="chartreuse3", "Mp"="#0099FF"),
                       labels = c("conservative (Mc)", "middle (Mm)", "progressive (Mp)")) +
    guides(color=guide_legend(title="Moral Position")) + 
    theme(legend.position="bottom",
          legend.text = element_text(size=10),
          plot.title = element_text(hjust = 0.5))
  
    data2 <- read_dta(data_name) %>% filter(`_predict`==3) %>%
    select('_at1', '_margin', '_ci_lb', '_ci_ub')
  E <- c("Ec","Em","Ep")
  data2$eco <- rep(E, 3)
  M <- c("Mc","Mm","Mp")
  data2$moral <- rep(M, each=3)
  colnames(data2) <- c("eco_moral", "pred", "ll", "ul", "eco", "moral")
  
  right <- ggplot(data2, aes(x=eco, color = moral)) + 
    geom_vline(aes(xintercept="Ec"), linewidth=30, color = alpha("gray", 0.2)) +
    geom_vline(aes(xintercept="Em"), linewidth=30, color = alpha("gray", 0.2)) +
    geom_vline(aes(xintercept="Ep"), linewidth=30, color = alpha("gray", 0.2)) +
    geom_pointrange(aes(y=pred, ymin=ll, ymax=ul), position=position_dodge(0.5)) +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"),
          axis.text = element_text(size = 12),
          axis.title.y = element_text(size = 14),
          axis.title.x = element_text(size = 14),
          axis.text.x = element_text(color = "black"),
          axis.text.y = element_text(color = "black")) +
    geom_hline(yintercept=overall[i,2], linetype="dashed", color = "red") +
    annotate("text", x=2, y=overall[i,2]+0.02, label=paste0("overall probability=", round(overall[i,2],2)), size=3, color="red") +
    ylab("Predicted Probability of Being Rightist") +
    xlab("Economic Position") +
    ggtitle("Right") +
    scale_x_discrete(breaks=c("Ec", "Em", "Ep"),
                     labels=c("conservative (Ec)", "middle (Em)", "progressive (Ep)")) +
    scale_y_continuous(limits = c(-0.1,0.9), expand = c(0, 0)) +
    scale_color_manual(values = c("Mc"="#FF99FF", "Mm"="chartreuse3", "Mp"="#0099FF"),
                       labels = c("conservative (Mc)", "middle (Mm)", "progressive (Mp)")) +
    guides(color=guide_legend(title="Moral Position")) + 
    theme(legend.position="bottom",
          legend.text = element_text(size=10),
          plot.title = element_text(hjust = 0.5))
  
  data <- rbind(data1, data2) %>% select(-eco_moral)
  data$ID1 <- rep(c("left", "right"), each=9)
  write.csv(data, file = paste0("model", i, "_pred_prob.csv"))
##########################################


  plot1 <- ggarrange(left, right, common.legend = TRUE, heights = 1,
                      legend = "bottom", ncol=2) + bgcolor("white")  
  ggsave(filename = paste0("model", i, ".png"), plot1, width=10, height=5, dpi=600)

}


